Анализ данных RRBS-seq

Введение

Метилирование ДНК играет важную роль в регуляции экспрессии генов млекопитающих. Метилирование промотора, как правило, инактивирует экспрессию группы рядом расположенных генов, так как транскрипционные факторы не могут связаться с промотором вследствие как прямого стереохимического влияния метильной группы, так и изменения пространственной структуры хроматина. ДНК метилируется по остаткам цитозина по атому углерода в 5 положении с помощью специальных ферментов ДНК-метилтрансфераз. Цитозины метилируются в большинстве своём в специальных областях, богатыми CG-динуклеотидами. Такие области называют CpG-островки.

Исследование метилирования генов в онкологии помогает понять природу изменений в опухолевой ткани, выделить подтипы опухолей в соответствии с которыми можно подбирать терапевтические агенты. Методов, измеряющих метилирование ДНК, изобретено большое количество. В этом мастер-классе мы поговорим о бисульфитном секвенировании с уменьшенной представленностью (Reduced Representation Bisulfite Sequencing, RRBS).

Внимание: методы и подходы, рассказываемые в этом блокноте также применимы к полногеномному бисульфитному секвенированию (Whole Genome Bisulfite Sequencing, WGBS). Данный блокнот не претендует на истину в последней инстанции. Существуют альтернативные подходы к анализу и соответствующие программы и библиотеки.

RRBS

RRBS - бисульфитное секвенирование с уменьшенной представленностью. Оно более дешевое, чем секвенирование полного генома, но даёт информацию более, чем о 80% CpG островках человеческого генома. Достигается этот эффект за счёт обогащения библиотеки специальными рестриктазами (MspI, XmaI).

Библиотеки в процессе подготовки к секвенированию подвергают бисульфитной конверсии, т.е. обрабатывают раствором NaHSO₃ . Суть этого шага, представлена на рисунке @ref(fig:bisulfite). Неметилированные цитозины при обработке бисульфитом натрия окисляются до урацила, который при ПЦР заменяется на тимин. Однако такому не подвержены метилированные цитозины. Следовательно, на том месте, где секвенатором получен тимин, был неметилированный цитозин, а процитанные цитозины - метилированные. Таким образом, можно измерять уровень метилирования в конкретном локусе ДНК.

Бисульфитная конверсия

Уровень метилирования рассчитывается по следующей формуле:

\[ methylation = \frac{C}{C+T} \]

Анализ данных

Последовательность обработки данных представлена рисунке @ref{fig:workflow} (Wreczycka et al. 2017).

Конвейер анализа данных RRBS

Контроль качества осуществляется с помощью стандартной программы FastQC.

library(DiagrammeR)
mermaid("
graph LR
    A{}-->B
")

Подготовка к работе

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install(c("GEOquery","limma", "edgeR", "methylKit", "DMRcate"))
install.packages(c("stringr", "plyr"))

Импорт пакетов

Загрузка данных

# Чтобы загрузка не падала из-за временного ограничения на обмен между сервером и клиентом (вами)
options(timeout = max(300, getOption("timeout")))
options(download.file.method.GEOquery = "wget")

experiment = GEOquery::getGEO("GSE122799")[[1]]

untar("GSE122799/GSE122799_RAW.tar", exdir = "GSE122799/")
files = Sys.glob("GSE122799/*.bismark.cov.gz")

sample.include = which(phenotypes$characteristics_ch1.2 == "mtypehe: hiHER" | phenotypes$characteristics_ch1.2 == "mtypehe: hiLumAB")
phenotypes1 = phenotypes[sample.include, ]
files1 = files[sample.include]

Описательная статистика

getMethylationStats(raw_data[[2]],plot=TRUE,both.strands=FALSE)

getCoverageStats(raw_data[[2]],plot=TRUE,both.strands=FALSE)

for (sample in raw_data) {
  getMethylationStats(sample,plot=TRUE,both.strands=FALSE)
  Sys.sleep(5)
}

filtered = filterByCoverage(raw_data,lo.count=10,lo.perc=NULL,
                                hi.count=NULL,hi.perc=99.5)

Объединение результатов

meth=methylKit::unite(filtered, destrand=FALSE,min.per.group = 10L, mc.cores = 30)
#load("Unite_meth.RData")

Фильтрация

methMatrix = percMethylation(meth,rowids = T)
sds=matrixStats::rowSds(methMatrix,na.rm = T)
print(dim(methMatrix))
hist(sds,col="cornflowerblue",xlab="Std. dev. per CpG", breaks = 50)

dim(meth)

meth.filt = meth[sds > 10,]

dim(meth.filt)

shared = with(
  list(xx = ifelse(is.na(methMatrix), 0, 1), n = ncol(methMatrix)),
  ldply(
    1:(n-1),
    function(i)
    {
      # построим список каждый-с-каждым
      ddply(
        data.frame("i"=rep(i, n-i), "j"=(i+1):n),
        c("i","j"),
        function(row){
          i <- row[1,1]; j <- row[1,2];
          sum(xx[,i] * xx[,j]);
        });
    },.parallel=TRUE))
print(shared[shared$V1 <= 10000,])

miss =  function(x){
  sum(is.na(x))/(sum(!is.na(x))+sum(is.na(x)));
}

meth.filt.mat = percMethylation(meth.filt)
miss_x <- miss(methMatrix);
max_miss = 0.2
while (miss_x > max_miss) {
  xx = ifelse(is.na(meth.filt), 0, 1);
  
  # в скольких образцах есть данные по каждому локусу
  cvrg_c = apply(xx, 1, sum);
  cvrg_c_mean = mean(cvrg_c);
  cvrg_c_sd = sd(cvrg_c);
  filter_out_plan = nrow(meth.filt)*(miss_x - max_miss);
  
  cvrg_c_thrsh = cvrg_c[tail(head(order(cvrg_c), filter_out_plan) ,1)];
  cvrg_c_flt = cvrg_c <= cvrg_c_thrsh;
  
  meth.filt = meth.filt[!cvrg_c_flt, ];
  meth.filt.mat = percMethylation(meth.filt)
  miss_x = miss(meth.filt);
  
  cat(sprintf(
    "Filters out: #plan=%f, #fact=%d CpGs  ==>  Dim: %d x %d (samples x CpGs); miss(x) = %.4f \n",
    filter_out_plan, sum(cvrg_c_flt), dim(meth.filt.mat)[2], dim(meth.filt.mat)[1], miss_x ));
}

Анализ после фильтрации

png(width = 2400, height = 2400, filename = "correlation.png")
getCorrelation(meth.filt, plot = T)
dev.off()

Корреляция уровня метилирования между образцами

PCASamples(meth.filt)

PCASamples(meth.filt, screeplot = T)

tree = clusterSamples(meth.filt, dist = "correlation", method="ward.D2")
rect.hclust(tree, k = 2, border = "darkgreen")

DMR

meth.diff = calculateDiffMeth(meth.filt, overdispersion = "MN", 
                              adjust = "fdr", 
                              test = "F",
                              mc.cores = 30)

meth.diff$dif.met = ifelse(meth.diff$meth.diff <= -20 & meth.diff$qvalue <= 0.05,
                           "HYPO",
                           ifelse(meth.diff$meth.diff >= 20 & meth.diff$qvalue <= 0.05,
                                  "HYPER",
                                  "NO"))
methylKit::bedgraph(meth.diff,"meth_diff.bedgraph",col.name = "meth.diff")

signif.diff = getMethylDiff(meth.diff,qvalue = 0.05, difference = 20)

diff.df = do.call(cbind.data.frame, meth.diff@.Data)
colnames(diff.df) = c("chr", "start", "end", "strand", "pvalue", "qvalue", "meth.diff", "dif.met")
ggplot(data=diff.df, aes(x=meth.diff, y=-log10(pvalue), col=dif.met)) +
  geom_point() + 
  theme_minimal() +
  scale_color_manual(values=c("HYPO"="blue", "NO"="black", "HYPER"="green")) +
  geom_vline(xintercept=c(-20, 20), col="red") +
  geom_hline(yintercept=-log10(0.05), col="red")

rownames = paste(signif.diff$chr,signif.diff$start,signif.diff$end,sep = ".")
counts.dif.met = percMethylation(meth.filt,rowids = T)[rownames,]
ann_col_df = as.data.frame(phenotypes1$characteristics_ch1.2)
colnames(ann_col_df) = "Phenotype"
ann_col_df$Phenotype = as.factor(ann_col_df$Phenotype)
levels(ann_col_df$Phenotype) = c("hiHER","hiLumAB")
row.names(ann_col_df) = phenotypes1$geo_accession
ann_col_df

ann_colors = list(Phenotype = c("hiHER"="red", "hiLumAB"="orange"))

pheatmap::pheatmap(counts.dif.met,
                   color = colorRampPalette(c("white","cyan", "blue"))(100),
                   clustering_distance_rows = "manhattan",
                   clustering_distance_cols = "manhattan",
                   clustering_method = "ward.D2",
                   scale = "none",
                   annotation_col = ann_col_df,
                   annotation_colors = ann_colors,
                   treeheight_row = 100,
                   treeheight_col = 75)

require(genomation)
gene.obj=readTranscriptFeatures("refseq.hg19.bed.txt")
gene.parts = annotateWithGeneParts(as(signif.diff,"GRanges"),gene.obj)


cpg.obj=readFeatureFlank("cpgi.hg19.bed.txt",
                         feature.flank.name=c("CpGi","shores"))
diffCpGann=annotateWithFeatureFlank(as(signif.diff,"GRanges"),
                                    cpg.obj$CpGi,cpg.obj$shores,
                                    feature.name="CpGi",flank.name="shores")

library("org.Hs.eg.db")
promoters = gene.parts@dist.to.TSS[gene.parts@members[,1] == 1,]
cols <- c("ENTREZID","GENENAME","SYMBOL", "GENENAME")
annotation = AnnotationDbi::select(org.Hs.eg.db, keys=unique(promoters$feature.name), columns=cols, keytype="REFSEQ")
coord_promoters = cbind.data.frame(promoters, signif.diff[promoters$target.row,])
annotated_promoters = coord_promoters %>% dplyr::inner_join(annotation, by=c('feature.name'='REFSEQ'))

write.csv2(annotated_promoters, file = "promoters.csv")

Rgraphviz

genesOfInterest = annotated_promoters$meth.diff
names(genesOfInterest) = annotated_promoters$ENTREZID
genesOfInterest = sort(genesOfInterest,decreasing = T)
go <- clusterProfiler::enrichGO(gene=annotated_promoters$ENTREZID, 
             ont ="BP", 
             keyType = "ENTREZID", 
             pvalueCutoff = 0.1, 
             OrgDb = org.Hs.eg.db, 
             pAdjustMethod = "fdr",
             readable = T)

head(go)

dotplot(go)

Wreczycka, Katarzyna, Alexander Gosdschan, Dilmurat Yusuf, Björn Grüning, Yassen Assenov, and Altuna Akalin. 2017. “Strategies for Analyzing Bisulfite Sequencing Data.” Journal of Biotechnology 261 (November): 105–15. https://doi.org/10.1016/j.jbiotec.2017.08.007.